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1. CONSERVATIVE DISSIPATIVE DIFFERENCE SCHEMES 


1.1. I ntroduction: Why Conservative Dissipative Difference Schemes? 

In solving problems of gas dynamics it is often permitted to ignore the 
dissipative processes in the gas, that is, viscous friction and heat 
conduction. This simplification of the physical picture boils down, mathe- 
matically, to a degeneration of the partial differential equations from 
second-order conservation laws, the Navier-Stokes equations, to first-order 
conservation laws, the equations of ideal compressible flow (ICF). Even in 
this approximation there remains a bewildering variety of complicated flow 
problems. 

Particularly notorious are the problems involving such a strong 
compression of the gas that, in spite of the a priori assumption, dissipation 
sooner or later dominates the flow, at last in certain regions known as 
shocks. In a shock the flow quantities undergo a significant change over a 
distance typical of the dissipative Interaction, l.e. the molecular mean free 
path. 

It is, of course. Impossible to infer the structure of a shock from the 
equations of ICF. The concept of ICF traditionally is saved and extended by 
representing a shock as a flow discontinuity. 

The motion of such an idealized shock may be derived from an Integral 
version of the first-order conservation laws, expressing the particular 
conservation principle for a finite volume of fluid and a finite lapse of 
time. Mathematically, this leads to the theory of weak solutions of nonlinar 
conservation equations. The algebraic equations relating the propagation speed 
of a dlscontinutiy to its amplitude are called the jump equations. 

Weak solutions are not unique. Shock discontinuties in which friction 
produces kinetic energy instead of dissipating it, and heat flows from lower 
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to higher temperatures, are as eligible as their physically realizable 
counterparts. Clearly a selection criterion must be invoked. We shall accept a 
weak solution of the first-order conservation equations only if it is the 
limit solution, for vanishingly small dissipation, of the second-order 
conservation equations. For gasdynamics this is equivalent to the following 
requirement: the entropy of the gas, measure of the accumulated effect of 
dissipation, must not decrease in a shock. This is called the entropy 
inequality. 

Thus, the advantage of lowering the order of the flow equations is partly 
offset by the need to introduce extra equations and an extra Inequality. In 
consequence, analytic treatment of ICF problems is impossible in all but a few 
cases. The numerical treatment, however, can be entirely successful. The key 
to success is the combination of conservation and artificial dissipation. 

The idea behind artificial dissipation is that, since in ICF the effect 
of dissipation is ignored, it may as well be e xaggerated . By providing a 
difference scheme for ICF with sufficiently large dissipative terms it is 
possible to achieve that shocks, whenever these appear, posses a structure 
coarse enough to be resolved in the computational grid. 

The concept of conservation enters if we make the difference scheme 
consistent with the integral form, rather than the differential form, of the 
equations of ICF. The scheme is then said to be "conservative" or to have the 
"conservation form"; another often-heard description is "finite-volume 
scheme". It has been shown, in theory and in practice, that, in using a 
conservative scheme, the numerical stability of a solution containing an 
admissible shock automatically guarantees the correct motion of the shock. But 
in order to achieve numerical stability and to select the proper shocks, the 
scheme indeed has to be dissipative. 
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1.2. Hist ory 

The first example of a conservative dissipative difference scheme was the 
first-order-accurate Lax- Friedrichs (LF) scheme, presented and analyzed by Lax 
(1954). It is the least accurate of its kind and has no practical value today. 
More influential was its second-order-accurate successor, the Lax-Wendroff 
(LW) scheme (1960) which still is widely used by aerodynamicists in the two- 
step form derived by MacCormack (1969). 

The survivor among first-order schemes is Godunov's (1959) method based 
on upwind differencing. The decision which direction is upwind is made on the 
basis of the speeds of the finite-amplitude waves by which discrete fluid 
volumes interact. With this strategy, the solution of Riemann's initial-value 
problem - the shock-tube or diaphragm problem - becomes a building block of 
the difference scheme. At present we see a rapid expansion of the literature 
on "approximate Riemann solvers"; Roe (1980), Osher (1980), Harten and Lax 
(1981), Van Leer (1982), Colella (1982); these find their way in first-order 
as well as higher-order upwind schemes. A review of this subject was given by 
Harten, Lax and Van Leer (1983); I shall return to it later in the present 
lecture. 

It turns out that the history of conservative dissipative difference 
schemes is divided more or less by the decades. 

During the fifties, first-order-accurate methods were developed and 
tested. These methods are what now is called "robust" and, schematically 
speaking, had no other problems than their modest accuracy. Particularly 
annoying was the numerical diffusion of entropy. 

The sixties gave us the second-order-accurate schemes, meaning different 
kinds of trouble. Numerical oscillations and nonlinear instabilities near 
shock waves largely spoiled the pleasure of the reduction in artificial 


dissipation. 
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During the seventies the design problems of the sixties were gradually 
solved. For Instance, the prevention of numerical oscillations is now well 
understood. As this knowledge has not been fully absorbed by the IGF 
community, I have made it the subject of the second lecture. 

The eighties have started lirtth an increased interest in explicitly using 
the physics embedded in the IGF equations for improving the numerical methods. 
In particular, there is a strong emphasis on upwind differencing. Some modern 
ideas on how to compute one-dimensional flow are covered by the present 
lecture; how well these ideas carry over to multi-dimensional flow is 
investigated in the third lecture. 


1.3. A family of second-order finite-volume schemes 

I shall Illustrate the recent developments in computing one-dimensional 
IGF on the basis of a family of second-order-accurate difference schemes. The 
initial-value representation and algorithm structure are the same as in Van 
Leer (1980); the notation is also the same, except for the present use of a 
superscript to Indicate the time level. 

We start with schemes for the scalar linear convection equation 


q + a q = o, a = constant. 


( 1 ) 


assuming a piecewise linear initial-value distribution 


qCx,!*') = q? + 


^iq 


1 Ax ^’^"^1^’ ^i-i = 


( 2 ) 


as in Van Leer (1980), Fig. 4. The gradient is computed by averaging the 
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numerical gradient values in the neighborhood of x^: 

'S.q” = i(I-k) a. + i(1+k) A. iq” . (3) 

1 1 ~ 2 iT"? 

For the moment, the weight k Is assumed to be independent of q. The three 
possible values of k most frequently met are 

K = sgn a = s, 0 and -s , (4) 

leading to the LW scheme, Fromm’s (1967) scheme and Moretti’ s (1979) 

X -scheme, respectively. 

Tlie evolution in one time-step of the initial-value distribution 
discretized as above is now computed exactly, namely, by shifting the 
distribution over a distance aAt. Averaging over the finite volumes then 
yields the updated values q^ , from which the gradient in each volume can 
again be determined. 

The LW scheme Involves only the downwind value of the gradient of q^. 

Shifting the inital-value distribution brings upwind information to x^, so 

that the LW scheme ultimately is Independent of the sign of a. The LW scheme 

is a so-called central-difference scheme and involves q^ , , q?, and q” , to 

i-1 1 1+1 

find . It is stable for 

|a| <1, (5.1) 

where 

|a|= aAt/Ax (5.2) 

is the Courant-Frledrichs-Lewy (CFL) number, and produces a predominantly 
negative phase error (see Fromm (1967)). 

Fromm's scheme is upwind biased, involves q^_2s» ^i-s’ *^i ^i+s’ 
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has a zero phase error at a = s/2 for all frequencies. Its stability domain is 
given by (5.1). 

The X-scheme, involving only the upwind value of the gradient of q”, 

remains fully one-sided; the values 1i_2s* '^i-s ^i used to 

n+1 

compute q^ . 

Its stability condition is 

\o\ < 2 , ( 6 ) 

but in practice we have to put up with (5,1), in order not to complicate the 
preservation of monotonicity. Under (5.1) the phase error of the X-scheme is 
predominantly positive (see Fromm (1967)). 


1.4. From co n vectio n to IGF 

In order to convert the above difference schemes for a single linear 
convection equation into schemes for the hyperbolic system of the nonlinear 
conservation laws of IGF, we must regard equation (1) as one of the 
characteristic (diagonalized) equations, q as a characteristic state quantity 
and a as a characteristic speed; see Van Leer (1980). Since a now depends 
on q, so does k, and it suddenly appears impractical (although feasible) to 
use any other value than k = 0. This leads to the following observation; 

In upwind-biased schemes the discretized initial-value distribution is 
Independent of the direction in which the various physical signals propagate ; 
see Van Leer (1977). 

Restricting ourselves to Fromm's scheme, we follow the recipe from [24] 
for advancing in time (the approach by Van Leer (1979) is different). We 
first update the distribution in each volume over a half step in time. 
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without regard to the neighboring volumes, with a central-difference scheme. 
For the characteristic quantities (other state quantities may be used) this 
bolls down to 





(7.1) 




(7.2) 


with sufficient accuracy. 

Subsequently, time-centered values are obtained on the upwind side of 
each interface. For ICF, described by the system of conservation laws 


w + [f(w)] = 0, 

u X 


( 8 ) 


this means we must compute the complete states w^ and w^^ on the left and right 
side of each interface and solve the local Rlemann problem, exactly or 
approximately, to find out what the upwind components are. In formula: 


n+^ _ 

"•(l+i)-! “ 


n+4 


± i«iq” , 



^(w 


n+i 

(i-H)- ’ 


’^(i-H)+^ ’ 


(9.1) 


(9.2) 


(9.3) 


here $(w ,w ) is the flux resulting at an interface after resolution of the 
L K 

initial discontinuity and F is the numerical flux vector used to approximately 
Integrate the system (8) over one time step and one finite volume: 


n+1 


w. 




5,n+i 

-±±L 


- F 


,n+i’ 

i-i 


= 0 


At 


+ 


Ax 


( 10 ) 
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Note that twice in this second-order scheme there is an exchange of 
information between neighboring volumes. 

The first exchange, finite-differencing of at t’^, is purely numerical 
- initial-value interpolation - and involves no physics. We may at this stage 
use any complete set of state quantities; the choice of characteristic 
quantities is imperative only near boundaries and advisable near discontinui- 
ties. 

The second exchange, solving the Riemann problem defined by w^ | i ^ - and 
is purely physical: interaction through waves. Simplifying the 

solution of the Riemann problem, while desirable for numerical reasons, is a 
matter of simplifying the physics of the wave Interactions. 


1.5, Approximate Riemann solvers 


I shall discuss two distinct approximate solvers of Riemann’ s problem: 
the one by Roe (1980), based on linear waves, and the one by Van Leer (1982), 
based on material transport. 

Roe solves the Riemann problem defined by the initial values 


w 


w. 


w. 


R 


X < 0 , 

X > 0 , 


by solving it for a linearization of Eq. (10) : 


( 11 ) 


w^ + A (w^, w^) = 0 , 


( 12 ) 
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A 

where A (w^^, Wj^) is some average of 

A(w) = df(w)/dw (13) 

over the interval (wj^, w^) . 

Mathematically it is clear that 

A(Wl,Wj^) = A( -- -) 4- 0((Wj^ -Wj^) ) , (14) 

which is fine for weak waves. Roe demands that the solution also be exact if 
there is only one wave with a nonzero amplitude; this wave, however, maybe 

A 

arbitrarily strong. This can be achieved only by chosing A such that 

A(w^,Wr)*(Wr-w^) = f(Wj^)-f(w^) , (15) 

the discrete analogue of (13). If w^^ and w^^ can be connected by a single 
discontinuity with speed U, the jumps in w and f are connected by the jump 
equation 

U(Wj^-Wl) = f(w^)-f(w^). (16) 

Upon comparison of (15) and (16) it follows that w^-Wj^ is an eigenvector of A, 
and U the corresponding eigenvalue. 

Harten (1981) has shown for hyperbolic systems (8) admitting an entropy 
condition, that the mean value 

A(w^,Wr) = Q/^A(w^+9(Wj^_Wj^))de (17) 

not only satisfies (15) but also has a complete set of real eigenvectors and 
eigenvalues, making Eq. (12) hyperbolic under all circumstance. 
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Roe (1981) has indicated an algebraic procedure to construct matrices 
A(w^jWjj^) of the form A(w) where w is some average of and His arguments 
bear heavily on the specific homogenlety properties of the conservation laws 
involved. 

Once A has been found, it is split in its positive and negative parts A"*" 

A 

and A ; positive and negative refer to the eigenvalues of the matrix. 

Correspondingly we may split Af s f(w ) - f(w ) in 

K L 

(Af)'*’ = a'*’ (w^ - w^) , (18.1) 

(Af) = A (w^ - w^) , (18.2) 

the changes of the flux across the forward moving and backward moving waves, 

respectively. This is called flux-difference splitting. 

It is easy to see that the flux $(w ,w ) needed in Eq. (9.3) equals 

L K 


HWl,Wr) = f(Wj^) + (Af) 


(19.1) 


= f(Wj^) - (Af)'*' ; 


(19.2) 


see also Apl, Fig. 8. 

A word of caution: the jump condition (16) does not imply the entropy 

condition: Wj^ and Wj^ may be exchanged. Likewise, if 

A A 

A 

the scheme incorporating A can not recognize inadmissible shocks and may not 
destroy these. Some asymmetry must therefore be introduced in the dependence 

A 

of A on Wj^ and v!-^; see Harten and Lax (1981). 
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An alternative to (18) is to split the flux Itself: 

f(w) = f^(w) + f (w) ; (21) 

f^(w) and f (w) may be called the forward and the backward flux, respectively. 
In the review by Harten, Lax and Van Leer (1983) it is demonstrated that this 
kind of splitting follows from interpreting the conservation laws (8) as 
moments of the collisionless Boltzmann equation, using a numerically 
convenient particle-velocity distribution. 

Van Leer (1982) shows how to achieve a splitting like (21) for the flux 
in the Euler equations. Care must be taken to ensure that f^(w) and 

f~(w) are continuously differentiable functions of w, so that Eq. (8) can be 

A ^ __ 

consistently approximated. In contrast, A (w^,Wj^) and A 
flux- difference splitting (18) just have to be continuous. In any case, 
second-order terms with split fluxes or split flux-differences must be 
avoided; see Mulder and Van Leer (1983). 


1.6 Alternative ways of time-stepping 


If temporal accuracy is not needed, for instance, in marching toward a steady 
state, we retain from the scheme from Sec. 1.4 only the steps relating to the 
approximation of f^. Following Jameson et al, (1981) we may advance in time 
with an m-step Runge-Kutta algorithm: 


(0) n 

' SS TJ 

1 i 

(k) n 

w. =* w. 

X i 


(k) 

5L 

k! 


At 


^Vax, k=l,..,m. 
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1 , 


(m) n+1 

7 . = W. 

1 1 


( 22 ) 


Once a steady state Is achieved It is independent of At. The same steady 
state would be found with the two-step scheme from Sec. 1.4, in the limit 
of At 0. 

Unfortunately, the stability condition associated with upwind-biased 

spatial-differencing operators is rather restrictive, due to the sizable 

negative real part of the eigenvalues of such operators. This is illustrated 

in Figure 1 (taken from work done in collaboration with E. Turkel (1982)), 

for the combination of Fromm's space-differencing and fourth-order four-step 

fk> 

Ringe-Kutta time-stepping (a^ = 1, k = 1,2, 3, 4). The stability conditon 

after four steps is a meagre 


a < 1 .4 ; 


(23.1) 


(3) (4) 

changing the coefficients a and a improves this up to 


I < 1*9 . 


(23.2) 


In contrast, the central-differencing operator of the LW scheme, which has 
purely imaginary eigenvalues, leads to a stability condition 

1 0 1 < 2.8 , (23.3) 

see Jameson et al. (1981). We conclude that upwind differencing is not 
efficient in combination with Runge-Kutta time-stepping. 

Another possibility is to use implicit time stepping, such as the 
"backward Euler" or "Implicit Euler" method. In practice the implicit 
difference scheme must be linearized in time to make inversion possible. The 
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linearization is only feasible if the numerical flux-function F is differen- 
tiable with respect to its arguments. This requirement favors the use of flux- 
vector splitting; see Mulder and Van Leer (1983). 


2. THE RECOGNITION AND REPRESENTATION OF DISCONTINUITIES 
2.1. Monotone init ial-value inte rpolation 

The formula (3) for is based on the assumption that the gradient of 

w varies smoothly between x^_^ and • If w or its first derivative is 

discontinuous in this Interval, the value of 6 thus computed is meaningless 
and, when used, may give rise to numerical oscillations. The sight of a shock 
with numerical oscillations is well known and needs no illustration here. 

A local discontinuity in w or w^ may be recognized, in a sequence of 
mesh-refinements, by a local extremum of ^.jiw/A ,w that goes to infinity or 
to a finite value 4 1. The only place in a smooth solution where this 

behaviour is also seen is at an extremum. 

It has been demonstrated by Van Leer (1972, 1974, 1977) that meaningful 
values of <5 are obtained under the extra constraint of monotonicity. For a 

linear convection equation it says that the updated discrete distribution must 
be monotone if the discrete initial-value distribution is monotone; this of 
course,, is true for the exact solution. Satisfying this condition turns out to 

simply be a matter of monotonlcly Interpolating the initial values inside 

VO lume s . 

To fix our thoughts, let 9^ ^ *^1+1 ^ ^ ^ shall use the 

abbreviations A = A. iq*^, A = A .q’’, d.q’^ = 6. It is easily seen that, for 
“* i”2 ‘ i *‘2 1 

n+1 

the preservation of monotonicity at t , we must limit 6 as follows: 
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6 < { 


T A, (no overshoot in x. , , for a ~ 0) 

1-0 + 1+1 

2 

— A (no overshoot in x. for a ~ 1) 
a - i 


(34) 


Note that, regardless of the value of o, 6 must vanish if A^ or A_ vanishes. 
Taking the minimum of condition (34) over all values of o, we arrive at 


6 < min (2A , 2A ) . 

— "T 


(35) 


With 6 within these bounds, the Initial values (2) inside volume i will not go 
beyond the neighboring zone averages; see Van Leer (1980). 

For the three schemes of Sec. 1,3, (35) means that no limiting is 

necessary if 


(LW) 0 < A_^/A_ < 2 , 


(36.1) 


(Fromm) 3 , 


(36.2) 


(X) i < A_^/A_ < ~ . 


(36.3) 


Using 


■l 


(A_^-A_)/(A^+A_) ~ ^ (In w^)^ 


(37) 


as a "smoothness monitor" (the logarithmic second derivative is a measure of 
curvature) , we can write (36) as 


(LW) -1 < < y , 


(38.1) 


(Fromm) “i ^ | , 


(38.2) 
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(X) - I < < 1 . (38.3) 

For the Fromm scheme, the permitted range (36.2) or (38.2) remains the 
same when the sign of a is flipped; this is not true for the ranges associated 
with the other schemes. Hence, the initial-value dicretization for Fromm's 
scheme remains independent of the physics. 

2.2. Limiters, switches and artificial dissipation 

Eq. (35) can be satisfied by limiting 6: 


6^. = min ( 2A , 6, 2A ) , 

11m - + ’ 


(39) 


but a smoother dependence on A /A is preferable. Harmful effects of clipping 

+ — 

caused by suddenly acting limiters like (39) are described by Van Leer (1977), 
Woodward and Colella (1982) and Van Albada et al. (1982). 

In the latter paper a smooth limiter for Fromm's scheme was described: 



(40.1) 

(40.2) 

(40.3) 


2 

here e is a small threshold that also prevents zerodlvide. The treshold can 


be chosen such that the limiting effect disappears in the neighborhood of a 

2 

smooth extremum, that is, when A^ and A have values ~ (Ax) . 
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Clearly, ~ (Ax)^ will do. 

Eq. (40.1) shows that is a weighted average of and A , biased 

toward the smaller of the two. In the limit of A^/A = “ or 0 , the average 
equals the smaller value. Eqs. (40.2) and (40.3) show that equals its 

unlimited value (3) multiplied by a switch factor; this factor falls about 

S 2 ~ KAx)^ (In (40.4) 

short of unity in a smooth solution. The appearance of the logarithmic second 
derivative is no surprise. 

A tighter fit to (39) is 


6 


lim 


(2A^ + e^)A_ + (2A^ + 

A? + 2|a,A + A^ + 2e^ 

+ I + — — 


(41) 


which tends toward twice the smaller of A_ and A^ , in the limit 

of A^/A = ” or 0. 

Eqs. (40) and (41) adso limit the value of S if the signs 

of A^ and A_ are different, that is, if the initial-value sequence is not 
monotone. The numerical experiments of Mulder and Van Leer (1983) suggest that 

(40) gives a cleaner representation of the solution near a steady shock than 

(41) . 

The effect of the switch factor in (40) on the scheme can be interpreted 
in several ways. According to (40.2) the switch turns off the nominal gradient 
in a zone; since it is the gradient that leads to second-order accuracy, the 
switch can be said to turn off the second-order term of the scheme. Symbolicly 
(ignoring the conservation form) : 
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^2 ^2 ^ 2 ^’ (^ 2 . 1 ) 

here L 2 and Lj^ are the second-order operator and embedded first-order 
operator. Since b]^~L 2 is predominantly a dissipative term, S 2 must be an 
artificial-dissipation coefficient. 

Writing (42.1) as 


(Lo) 


2"^ lim 


(1 - S„) 




S2 Lj; 


(42.2) 


reveals that (^ 2 )]^^^ symbolizes a so-called hybrid scheme, see Harten (1978). 
Writing (42.1) as 


(L2)iim = + (1 - S2)(L2 - L^) (42.3) 

reminds us of the Flux-Corrected Transport (FCT) methods of Boris and Book 
(1973), where a first-order result is improved by adding a limited second- 
order term. These, however, are the methods mentioned earlier that show the 
disastrous clipping effect. 


2.3 Nonlinear case 


The analysis in 2.1 shows that for a linear convection equation the 
monotonicity of the numerical solution can be preserved by monotone 
interpolation of the solution inside each volume. This recipe works 
surprisingly well for nonlinear convection and seldom needs adjustments; see, 
however. Woodward and Colella (1982), Colella and Woodward (1982). For single 
nonlinear conservation laws it can be made exact; see Harten (1982). 
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What do shock profiles, computed with monotonicity preserving schemes, 
look like? To answer this question an example of shock propagation by three 
upwind-biased schemes based on monotone initial-value interpolation is given 
in Figure 2, taken from Woodward (1980). 


2 .4 Related subjects 


So far we have only discussed the role of the initial-value 
representation in preserving raonotonicity. Obviously the choice of the time- 
stepping algorithm will Influence the numerical results. The sharpness of 
steady shocks, for Instance, depends on the Riemann solver used, although the 
differences among second-order schemes are much less important than among 
first-order schemes; see e.g. Van Leer (1981)., (This subject will be reviewed 
in a forthcoming paper by Enqulst, Osher, Roe, and Van Leer (1983)). 

Discontinuities that have spread too much, in particular contact 
discontinuities, may be steepened by applying Harten’s (1977) artificial 
compression method (ACM). 

It is possible to use a smoothness monitor for other purposes than to 
limit gradient values. Some applications are: to indicate where the grid 
should be refined or moved, and where shock-fitting or shock-recovery should 
be carried out; for the latter subject see Morton (1982). 

3. MULTI-DIMENSIONAL METHODS 
3 . 1 Multi-dimensional convection 

In order to illustrate how the one-dimensional schemes of Section 2 can 
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be applied in multi-dimensional problems, I shall discuss two first-order 
schemes for convection in two dimensions. 

Godunov's one-dimensional scheme is what we get if, in the schemes 
studied previously, 6q is always set equal to zero. In other words, the 
initial values in each volume are chosen to be uniform, which accounts for the 
loss of the second-order accuracy. The Lax-Friedrichs scheme also uses uniform 
initial values in each volume but the grid is staggered in time, I shall start 
with the latter. 

The convection equation to be approximated is 

qt + a>0,b>0. (43) 

For convenience the points (x^, y^), y^), (xj,, Yj_ 2 ) and Yj-^) > 

and the volumes centered at these points are indicated by A, B, C and D, 
respectively; the center (xj^_i, Yj_i) of staggered volume is indicated by 

M. This is shown in Figure 3a. 

After moving the piecewlse-unlf orm distribution with a velocity (a,b) 
over a time interval At we find, by weighting q according to area (see Figure 
3b): 


^Tur'“ ”■ “CJ )Q? + 

^ M X y ^A X y 

+ (2 ~ + ®y)qc + ^2 + '^x’^ ^ ^ 


For comparison with Eq. (10) we imiy rewrite this as 


n+1 1 / n , n , n , tk 

I M ■ + "b + “c + ■’d’ 
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n+7 


V’ E 



) - a^(, 


n+| 

N 



(44.2) 


Here q is the time-average of the space average of q on the East side (AC) 

of volume M; the other averages are defined similarly. This formula differs 

from the version given orignally by Lax (1954), where the fluxes through the 

sides of volume M are not centered in time. This removes the cross-difference 

term oa(q^-q^-q^ + q^) from the scheme and reduces the scheme's 
X y ^A D 

stability domain. 

Scheme (44) can be written as the product of two one-dimensional LF 
schemes: 


- [Ki+Ty-h - yi-Ty-bitia+T/') - »,a-T/biq". 

E L^^(At) L^^(At)qJj , (45) 

where and denote a forward translation by one volume in the x- and y- 
direction. Hence, the method of fractional time-steps or time-splitting is 
exact for the linear LF scheme. 

Godunov' s scheme, illustrated in Figure 4, can be written as 


nd" 1 f ^ \/i 

q 4 - (i-q^)a-Oy)q 


A 


+ 0 ( 1-0 )q" + o ( 1-0 )q" + o o q" (46.1) 
X y B y X C x y D 


or 
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a [{(Ma^)q"44a^q"} - {(.Wo^)qc + 



/ n+i 
0 (q ^ 
X ^ E 



) - a^(q 


n+^ 

N 


n-H-. 
S ^ 


(46.2) 


where E, W, N and S now refer to the sides of volume A. Again the proper 
cross-difference tem is included and exact factorization is possible: 


n+l 


-1 


= {(1-0 (1-T 

^ y y 


)} {l-0^(l-T^b}q 


n 

ij 


E Ly(At) L^(At) 



(46) 


3.2 From convection to ICF 


The LF scheme can be exactly Implemented for a nonlinear system of 
conservation laws 


w^ + [f(w)]^ + [g(w)]y = 0 . 


(47) 


The 


flux aq 


nti 

E 


through the East side of M is changed into F 


n+4 

E ’ 


with 


== L _ 

E AtAy 



wJJcCt.y)) dy dt , 


(48) 


where w^^ (t,y) is the time-dependent solution of the one-dimensional Riemann 
problem on the line CA with initial values at t’'^. The other fluxes are 
obtained similarly. 

An useful approximation of (48) results if we first average w: 
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l_ 

At Ay 



w^^(t,y) dy dt 


= i(w? + 




At , n n. 
2Ax ’ 


( 49 ) 


and then adopt 



f (w' 


n+2 ^ 


9 


(50) 


which differs 0((Ax)^) from (48). Note that, owing to conservation, the 

ndv 

computation of the average w does not require knowledge of the detailed 

E 

Riemann solution. 

Eq. (49) is, of course, the one-dimensional LF scheme for flow along the 
line CA. This inspires the following Implementation of the LF scheme for Eq. 
(47). 


~n+t 


^LF^At,. n 

S '“2) ”13 


(51.1) 


~nfi 

w. 

i-il 


^LF,At, n 
”13 


(51.2) 


n+1 1 . n , n , n , n 

i-tl-i 4^ ij i~lj ij -1 i-lj-r 


At 

Ax ^ ij-| 


^i-lj-i^ Ay ^®1-K1 




(51.3) 


This version is symmetric in x and y, just as the time-splitting suggested by 
Strang (1968) for two-dimensional operators: 


L2p(At) = |{L,x(^t)L„(At) + L,/At)L^(At)} , 


X 


(52.1) 
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It must be said, however, that computational fluid dynamists almost 

exclusively use pure product operators like (45). To reduce the asymmetry 

betwen x and y, the sequence L L is alternated with L L : 

y A. y ^ 


L^(2At) = L (At)L (At)L (At)L (At) . 
2 X y y X 


(52.2) 


The exact formulation of Godunov's scheme for Eq. (47) requires knowledge 
of the solution of Riemann's problem in two dimensions, with initial values 
uniform in each quadrant. This solution is only known for a single 
conservation law and was obtained by Wagner (1980). In the same way as (51) is 
derived from (44b), we may derive from (46.2) an approximation to Godunov's 
scheme based entirely on the solution of one-dimensional Riemann problems: 


~n+i 

w. . 

id 


G At. n 

K (^’"u . 


(53.1) 


~nf| _G /At. n 
w. . = L (— r)w. . , 

ij X 2 11 ’ 


(53.2) 


w 


n+1 

Ij 


= w. 


11 


Ax 




_ 9^+^ 1 


4t ^;:rn+i 
Ay ^^ij+i 


- 

IJl 2 


(53.3) 


Here 


^+2 

^i+ii 


$(w 


•n+i 
11 


~n+^ . ;?n+j _ i/^n+f 


w 


;n+^ , 
i 1+1 ' 


(53.4) 


with the Blemann flux i>(w ,w ) defined previously. Thus, the one-dimensional 

L R 

Riemann solver is used twice in each space dimension. 

Q 

Expressing (53) in terms of only leads to 
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L2o(At) = I + (L^(At) - I) Ly(^) 

+ (Ly(At) - I) L^(^) . (54.1) 

Q 

Since Ljp(x) is linear in t, this may also be written as 

L2p(At) = I + (L^(At) - I) i(Ly(At) + I) 

(L^(At) - I) i(L^(At) + I) . (54.2) 

y X 

Using the homogeneity property of the Euler equations and, therefore, of the 
Rlemann solution, we may further work this out: 

L2p(At) = |[L^(At) (L^(At) + I) - L^(At). 

+ Ly(At) (L^(At) + I) - Ly(At)] , (54.3) 

Q 

This would reduce to Strang’s splitting (52.1) if would have the 
distributive property. This, of course, is only true in the linear case. 

3.3. Second-order schemes 


The examples of the previous section teach us that it is possible to 
mimic two-dimensional wave-interactions by combinations of one-dimensional 
wave interactions. These schemes contain terms ~ (At)^, but only have first- 
order accuracy. It is possible to apply the ideas of Section 3.1 and 3.2 to 
the second-order schemes of Section 2. This clearly Involves more algebra and 
will not be carried out here. One of the reasons is that none of these two- 
dimensional schemes has yet been programmed. It remains to be shown in 
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practice if the more subtle ways of splitting like (53) are indeed preferable 
to a pure product operator like (52.2). 


3.4. Related subjects 

In two dimensions the grid plays an even more Important role than in one 
dimension. Various ingredients of the schemes discussed in these lectures may 
be used more than once, in order to make decisions about refining the grid, 
moving the grid or just rotating the frame in which the solution of a one- 
dimensional Riemann problem is considered. 

For the latter decision we may use the direction of the gradient of the 

solution - readily available from 6 w and 6 w in a second-order scheme - or 

X y 

the projection angle at which the initial values in a Riemann problem could be 
best connected ("best" in some least-squares sense) by a single running wave. 
The latter idea is due to Harten (private communication). 

Moving a grid interface at the speed of the best-fitting single wave was 
realized in one dimension by Harten and Hyman (1982); two-dimensional results 
will soon follow. It appears that the Improvement due to the grid movement is 
dramatic for a first-order scheme and not nearly so much for second-order 
schemes. 

A two-dimensional extension of Roe's matrix theory is being considered 
(Baines, M. J., (1982), University of Reading, United Kingdom, private 
communication) . 
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Fig. 1. Stability domain of four-step Runge-Kutta methods in the complex 
plane. 

, (1) (2) (3n (4-v 

a,) With a =a =a = a =1. The inner contour is the locus of the 


Fourier transform of the second-order upwind spatial-differencing operator, 


for 0 = 1.3. 

(1) (2) (3) (4n 

b.) With (x^ = a = 1, a = .86, a = .44; a = 1.9. 


VELOCITY 


Fig. 2. Propagation of a one-dimensional 
shock, as represented by first-order (a), 
second-order (b) and third-order (c) up- 
wind methods. Note that the overall dis- 
hribution is always monotone. Reproduced 
from Woodward (1980). 










Fig. 3. Two-dimensional convection by the Lax-Friedrichs scheme, (a) Staggered 
grid of finite volumes, (b) Contours of the convected distribution after one 
time step (broken lines). The updated average value in zone M is an area- 
weighted average of the values in A, B, C and D; the areas are shaded 
distinctively. 
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